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Abstract 

In this project 1 , we study the hidden Markov ran- 
dom field (HMRF) model and its expectation-maximization 
(EM) algorithm. We implement a MAT LAB toolbox named 
HMRF-EM-image for 2D image segmentation using the 
HMRF -EM framework 2 . This toolbox also implements 
edge-prior-preserving image segmentation, and can be eas- 
ily reconfigured for other problems, such as 3D image seg- 
mentation. 

1. Introduction 

Markov random fields (MRFs) have been widely used 
for computer vision problems, such as image segmentation 
[4], surface reconstruction [ ] and depth inference [2]. 

The HMRF-EM framework was first proposed for seg- 
mentation of brain MR images [ ] . Given an image y = 
(2/1, . . . , yjq) where each yi is the intensity of a pixel, we 
want to infer a configuration of labels x = (x±, . . . , xn) 
where xi £ L and L is the set of all possible labels. In a bi- 
nary segmentation problem, L = {0, 1}. According to the 
MAP criterion, we seek the labeling x* which satisfies 

x* = argmax {P(y|x, 0)P(x)}. (1) 

X 

The prior probability P(x) is a Gibbs distribution, and the 
joint likelihood probability is 

P(y|x,9) = Y[P( yi \x,G) 

i 

= HP(yi\x i: 8 Xi ), (2) 



This work originally appears as the final project of Prof. Birsen 
Yazici's course Detection and Estimation Theory at RPI. 

2 This toolbox can be downloaded at the author's homepage 
http : //homepages .rpi . edu/ wangq 1 . 



where P(yi\xi, 6 Xi ) is a Gaussian distribution with parame- 
ters 6 X . = (n Xi ,cr Xi ). = {6i\l £ L} is the parameter set, 
which is obtained by the EM algorithm. 

2. EM Algorithm 

We use the EM algorithm to estimate the parameter set 
6 = {0i\l £ L}. We describe the EM algorithm by the 
following: 

1. Start: Assume we have an initial parameter set 6^. 

2. E-step: At the tth iteration, we have 6^, and we cal- 
culate the conditional expectation: 

Q(0|0 (£) ) = p[lnP(x,y|0)|y,0^~ 

= ^P(x|y,0«)lnP(x,y|0),(3) 

where x is me set of a U possible configurations of la- 
bels. 

3. M-step: Now maximize Q(Q\Q^) to obtain the next 
estimate: 

9 (£+1) = argmax Q(6|6 W ). (4) 
e 

Then let 0( £+1 ) — >• 0W and repeat from the E-step. 

Let G(z; 61) denote a Gaussian distribution function with 
parameters Q\ = (/i/, 07): 

We assume that the prior probability can be written as 

F(x) = |exp(-C/(x)), (6) 
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where U (x) is the prior energy function. We also assume 
that 

^(y|x,e) = Y[p(yi\xi,0 Xi ) 

i 

= Y[G( yi ;9 Xi ) 

i 

= ^exp(-t/(y|x)). (7) 

With these assumptions, the HMRF-EM algorithm is given 
below: 

1 . Start with initial parameter set 0(°) . 

2. Calculate the likelihood distribution P^ (yi\xi,9 Xi ). 

3. Using current parameter set 6^ to estimate the labels 
by MAP estimation: 



At) = 



argmax {P(y|x, 6 (£) )P(x)} 



= argmin {U(y\x, 6«) + U(x)}. (8) 

The algorithm for the MAP estimation is discussed in 
Section 3. 

4. Calculate the posterior distribution for all I G L and all 
pixels yi : 



P (t) (l\Vi) 



(9) 



where a;^. is the neighborhood configuration of xf', 
and 
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Note here we have 



p(i\4l) 




V c {l,xf) I .(11) 



5. Use (l\yi) to update the parameters: 

T.P (t \i\yi)vi 



(*+i) 



Mi 



(*+l)\2 



(12) 



E^lifcXifc-/^) 2 



-■(13) 



3. MAP Estimation 

In the EM algorithm, we need to solve for x* that mini- 
mizes the total posterior energy 



x* = argmin {J7(y|x, 9) + U(x)} 

with given y and 0, where the likelihood energy is 

C/(y|x,9) = Y^U( yi \xi,Q) 

i 

\2 



(14) 



In a T 



(15) 



The prior energy function U (x) has the form 

Z7(x) = ^K(x), (16) 

cec 

where V c (x) is the clique potential and C is the set of all 
possible cliques. 

In the image domain, we assume that one pixel has at 
most 4 neighbors: the pixels in its 4-neighborhood. Then 
the clique potential is defined on pairs of neighboring pix- 
els: 



V c (Xi, Xj) — ^ (1 Ixi,Xj)i 



where 



if X ^ X j 

if X> o X/ n 



(17) 



(18) 



We have developed an iterative algorithm to solve (14): 

1. To start with, we have an initial estimate x^ ^ , which is 
from the previous loop of the EM algorithm. 

2. Provided x^ k \ for all 1< i < N, we find 



„(*+!) 



argmin {Ufa \l) + ^ V c {l,xf ] )}. (19) 
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3. Repeat step 2 until J7(y|x, 6) + U(x) converges or a 
maximum k is achieved. 

4. Edge-Prior-Preserving Image Segmentation 

To use HMRF-EM framework for image segmentation, 
first we generate an initial segmentation using k-means 
clustering on the gray-level intensities of pixels. The initial 
segmentation provides the initial labels x^ ^ for the MAP 
algorithm, and the initial parameters for the EM al- 
gorithm. Then we run the EM algorithm, and the resulting 
label configuration x will be a refined segmentation result. 

Now we would like our segmentation to preserve the 
edges obtained by some edge detection algorithm. Assume 
we have a binary edge map z, where Z{ = 1 if the ith pixel 
is on an edge, and Z{ = if not. Then we modify (19) to 

xf +1) = argmin {Ufa\l) + ^ V^xf)}. (20) 

leL jeNi, Zj =o 
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Figure 1: Edge-prior-preserving image segmentation results, (a) Original image, (b) Canny edges, (c) Gaussian blurred 
image, (d) Initial labels obtained by k-means, where k = 2. (e) Final labels obtained by HMRF-EM algorithm, (f) Total 
posterior energy in each iteration of the EM algorithm. 



5. Experiment Results 

We run our HMRF-EM edge-prior-preserving segmenta- 
tion algorithm on example images. The binary edge map z 



is obtained by performing Canny edge detection [1] on the 
original image, and the observation y is obtained by per- 
forming Gaussian blur on the original image. Some results 
are shown in Figure 1. We can see that the initial labels 
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File 


Type 


Usage 


demo . m 


Runnable script 


A demo showing how to use the toolbox. 
Users can run this file directly. 


image _kme an s . m 


Function 


The k-means algorithm for 2D images. 
This will generate an initial segmentation. 


HMRF_EM . m 


Function 


The HMRF-EM algorithm. 


MRF MAP . m 


Function 


The MAP algorithm. 


U_YX . m 


Function 


Calculating the likelihood energy. 


U_X.m 


Function 


Calculating the prior energy. 


U_l .m 


Function 


Calculating the clique potential. 


gauss ianBlur . m 


Function 


Blurring an image using Gaussian kernel. 


gauss ianMask . m 


Function 


Obtaining the mask of Gaussian kernel. 


ind2i j .m 


Function 


Index to 2D image coordinates conversion. 


G.m 


Function 


Calculating Gaussian distribution probability. 


BoundMirrorExpand .m 


Function 


Expanding an image. 


BoundMirrorShrink .m 


Function 


Shrinking an image. 


Beijing World Park 8. JPG 


Image 


An example input image. 



Table 1: Name and usage of each file in the HMRF-EM- image toolbox. 



obtained by the k-means algorithm are not smooth enough, 
have morphological holes, and do not preserve the Canny 
edges. The HMRF refined labels overcome all these disad- 
vantages. 

6. Toolbox Documentation 

We provide the name and usage of each file in our 
MAT LAB toolbox HMRF -EM- image in Tabel 1. The 
U_X . m file can be modified to re-define pixel neighbor- 
hood relationships, and the U_l . m file can be modified to 
re-define the clique potentials. To reconfigure this toolbox 
for 3D image segmentation, the indexing system must be 
modified in several files. 

7. Discussion 

Our HMRF -EM- image toolbox is an implementation of 
the hidden Markov random field and its EM algorithm. This 
toolbox is well commented and easy to reconfigure. We 
have demonstrated the effectiveness of our toolbox on a 
simple example image. The HMRF model is mainly used 
to refine the direct segmentation output of some other algo- 
rithms. To get better segmentation results on more compli- 
cated images, some higher-level features should be used to 
construct y instead of just pixel intensities, and some more 
advanced algorithm should be used to generate the initial 
labels. 
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